Exploiting Phase Transitions in Catalysis: Adsorption of CO on doped VO2‐Polymorphs

Abstract VO2 is well known for its low‐temperature metal‐insulator transition between two phases with tetragonal rutile and monoclinic structure. The adsorption of CO on the two polymorphs of Mo‐doped VO2 is calculated to investigate the effect of a substrate phase change on the adsorption energy. The system is investigated theoretically at density‐functional theory level using a hybrid functional with London dispersion correction. We establish a computational protocol applicable for the study of physisorption on open‐shell transition metal oxides. The main task is to control the spin state of open‐shell slab models used to model adsorption of closed‐shell molecules in order to obtain numerically stable adsorption energies and to reduce spin contamination within the broken‐symmetry unrestricted Kohn‐Sham approximation. Applying this procedure, it is possible to identify the most stable adsorption positions of CO on both phases of VO2. CO adsorbs vertically with the C atom on a surface V atom in the monoclinic phase with an adsorption energy of −56 kJ/mol. The same adsorption position has an adsorption energy of only −46 kJ/mol on the rutile phase. Similar differences were obtained with multireference methods using an embedded cluster model. This effect may inspire experimental strategies exploiting the rutile↔ monoclinic VO2 phase transition in catalytic processes where CO is formed as product or as an intermediate.


Introduction
To meet the challenges of overall higher energy demand and the need for more sustainable resources faced by heterogeneous catalysis, many approaches to optimize industrial processes have been explored. [1,2] For example, fluctuating reaction conditions in catalytic process are investigated to improve reaction conditions. [2] In this way desorption of the product is facilitated, which would otherwise block the catalytically active sites. The present study investigates a phase transition of the substrate as an example of fluctuating reaction conditions. In order to improve the efficiency of a catalyst, the adsorption properties on the different phases of the catalyst should vary significantly. Then it would be possible to remove a product formed on the surface of one catalyst phase after transition to another phase which is less attractive.
In the present study we chose vanadium dioxide VO 2 as substrate. Vanadium oxides have been investigated as catalyst materials because of their rich redox chemistry. [3] VO 2 in particular is catalytically active for the desulfurization of dibenzothiophene, [4] the oxidative dehydrogenation of propane [5] or the electrochemical reduction of trinitrotoluene. [6] VO 2 shows a metal-insulator transition at 340 K between the tetragonal rutile and a monoclinic phase. [7,8] The effect of this phase transition has been investigated for the catalytic oxidative desulfurization of dibenzothiophene. [4] In the present study, we chose the adsorption of CO as a simple model system. CO is a product or an intermediate in many catalytic processes. [9,10] The adsorption of CO on the most stable surfaces of both VO 2 polymorphs is calculated in order to identify differences in the interaction strength which could be exploited in fluctuating reaction conditions for easier removal of the molecule.
Suitable density-functional theory (DFT) methods for modeling the energetic, electronic and structural properties of the insulating monoclinic (M 1 ) and the metallic rutile (R) phase of VO 2 were identified in previous studies. [11,12] A self-consistent hybrid functional (sc-PBE0) was found to provide the most accurate results for relative phase stability and electronic structure. Subsequently the surface properties of both VO 2 polymorphs were investigated. [13] A reconstruction of the R phase surfaces was found which has also been found experimentally. [14] Doping with Mo stabilizes the R phase and preserves the structure of both phases also at the surface. The most stable surface models with the Mo dopant in the outer layers of the surface models were used in the present study.
The adsorption of molecules on surfaces containing openshell transition metals has to be calculated with spin-polarized Kohn-Sham DFT. [15] This applies to the present system, since the VO 2 M 1 phase has an antiferromagnetic ground state according to Quantum-Monte Carlo calculations, [16] and the R phase is paramagnetic. [16][17][18] However, the mechanism of the metalinsulator transition and the electronic structures of the phases are not fully understood by theory and experiment, as shown in a review article by Poguet. [7] Electron-electron interactions were found to play a decisive role in the transition mechanism. This underlines the importance of spin polarization in the theoretical treatment of this system. In this study, we show the importance of controlling the spin state for a numerically stable calculation of adsorption properties. Additionally, the effect of spin contamination due to the broken-symmetry approximation is considered, which was discussed in previous studies. [19]

Computational Details
All calculations were performed with the crystal-orbital program CRYSTAL17 (version 1.0.2). [20] The hybrid functional M06 with D3 dispersion correction [21,22] was applied instead of the previously established sc-PBE0 method. [12,13] This was advantageous since M06 showed higher numerical stability in the adsorption calculations. Additionally, optimized D3 parameters are available for M06 which is not the case for sc-PBE0. Inclusion of London dispersion was expected to be crucial for the calculation of adsorption energies. Results for the bulk phases of VO 2 obtained with M06-D3 and sc-PBE0 are compared in the supplementary material Tables S1 and S2. M06-D3 is only slightly less accurate than sc-PBE0 for the calculation of lattice parameters, relative stability, and electronic structure.
Pob-TZVP-rev2 basis sets [23] were used for V, O and C, and a pob-TZVP basis set [24] for Mo. The present basis sets are larger than those of our previous study. [13] This proved to be necessary to ensure stable spin states and to reduce the basis set superposition error in the adsorption calculations. The integral truncation tolerances were set to the standard values of (10 À 6 , 10 À 6 , 10 À 6 , 10 À 6 , 10 À 12 ) to reduce computational effort. A Monkhorst-Pack net with 12 × 12 × 1 k-points was applied for the surface unit cells which was considered as converged. The optimized bulk lattice constants in Table S2 were used for the construction of surface models. To reduce computational cost only four-layer slab models of Modoped VO 2 were used. The most stable surfaces R (110) and M 1 (011) were investigated with one Mo atom each in the outermost layers of the slab. The R (110) surface is calculated with a 2 × 1 supercell which has the same size as the M 1 (011) primitive unit cell. Both surface unit cells contain four metal and eight oxygen atoms per layer. The optimized bulk lattice constants were applied for the construction of surface models. Surface models were relaxed by optimization of all atomic positions under symmetry constraints with fixed lattice vectors, applying the quasi-Newton algorithm implemented in CRYSTAL17 [20] unless indicated otherwise.
Ferromagnetic (FM), ferrimagnetic (FI) and antiferromagnetic (AFM) spin states were investigated for the surface models. Several methods were tested to set up the initial spin configuration of a system prior to the self-consistent cycle. The keyword SPINLOCK constrains the number of electrons with up and down spin (n α -n β ) until a predefined SCF convergence criterion (in the present calculations 10 À 5 a.u.) is reached. The spin states obtained in this way will be denoted as fixed spin(x) with x = n α -n β . An AFM state corresponds to x = 0. The keyword FDOCCUPY allows to define individual initial occupation numbers for the d orbitals of V and Mo. For VO 2 , it was found that an initial occupation of the d x 2 À y 2 orbital leads to the most stable spin states. These calculations were denoted as FDO.
A counterpoise correction was applied to the calculated adsorption energies E ads by replacing either the molecule or the slab atoms by ghost functions in single-point energy calculations of the optimized geometry. The correction is in most cases relatively small, � 10 kJ/ mol, due to the use of BSSE-corrected pob-TZVP-rev2 basis sets.
FI and AFM states can only be set up using the broken-symmetry approximation due to the use of single-determinant wavefunctions in DFT. Since this is a crude approximation of the exact wavefunction, we additionally performed complete active space SCF calculations for comparison. CASSCF and NEVPT2 calculations were performed for finite embedded clusters shown in Figure 1 with ORCA version 4.2.1. [25,26] Only the adsorption on V was considered in the comparison. The cluster calculations used def2-TZVPP as well as def2-QZVP [27] basis sets and the auxiliary basis def2/JK. [28] The clusters were embedded in a field of V-atom ECPs [29] and point charges. The point charge values were adjusted to the Mulliken charge of the central V atom, resulting in values + 1.8 and À 0.9. In the CASSCF calculations the five d-orbitals of the V-atom with one unpaired electron in the + IV oxidation state were considered to be in the active space. Only the lowest active CAS-orbitals were considered in the NEVPT2 calculation.

Results and Discussion
During preliminary calculations we found that the spin states of the slabs with and without adsorbate differed in an unpredictable way. Consequently, the calculated CO adsorption energies scattered significantly. Since the interaction of CO with V or Mo is weak to moderate, these changes in the spin state of the substrate are considered as artifacts of the self-consistent field procedure. Additionally, adsorption energies obtained from closed-shell calculations were found to vary significantly, even after small changes of the adsorbate structure (see Table S3). Therefore, we used various methods to control the spin configuration of the surface models and its effect on the adsorption energy (see below). The aim was to converge the slab to the same spin configuration with and without adsorbate, which should correspond to the magnetic ground state in both cases. The adsorption energies E ads are calculated for CO adsorbed C-down on a metal atom (see Figure 2). The resulting E ads and the relative energy ΔE = E spin À E groundstate to the spin state closest to the ground state are calculated and shown in Table 1.
We used various magnetic configurations of the Mo 2 V 14 O 32 slabs (see Table 1) with a fixed geometry (see Figure 2). In the formal + IV oxidation state, each V atom has a d 1 configuration and Mo has a d 2 configuration. Therefore the highest sensible  number of unpaired electrons is 18. We also defined an AFM spin configuration where the symmetry equivalent V and Mo atoms had opposite spin (fixed spin(0)). FI configurations with x = 8 and x = 16 were considered for comparison. The x = 16 state is the most stable spin state for undoped VO 2 surfaces. Some surface models were found to converge to a x = 8 spin state where the three V-atoms and the Mo-atom in the topmost layers showed an antiferromagnetic configuration while the inner layers showed ferromagnetic ordering. Therefore, the x = 8 state was also chosen. Manual spin setting for individual metal d orbitals (using the CRYSTAL keyword FDOCCUPY, denoted as FDO) was tested without spin locking, and with fixed spin(0) and fixed spin (18). The FM fixed spin (18) wavefunction with FDO was found to be the most stable spin configuration in both phases with and without the adsorbate. All other attempts resulted in SCF solutions that were up to 310 kJ/mol higher in energy. The stability of the spin states has a large influence on the adsorption energy. An unstable spin state in the surface without adsorbate leads to a large E ads . This can be seen in the rutile phase for the fixed spin(8) state. The surface is 310 kJ/mol less stable than the fixed spin(18) state with FDO, while the model with adsorbate was 116 kJ/mol less stable. After counterpoise correction this leads to a strongly negative E ads of À 93 kJ/mol. This can also be observed in the M 1 phase model calculated with FDO without fixed spin. The surface in this state is 265 kJ/mol less stable than the fixed spin(18) state with FDO. This leads to E ads of À 107 kJ/mol. On the other hand, an unstable spin state of the surface with adsorbate leads to a small or positive adsorption energy. This can be observed in the AFM fixed spin(0) and FDO (w/o fixed spin) state of the R phase. These states yield E ads of + 88 and + 14 kJ/mol, respectively. These results show the importance of finding the most stable spin state when calculating adsorption energies of open shell systems.
A previous study of Biz et al. [15] showed the importance of spin-polarization in calculations of surface adsorption due to exchange interaction. Our results show, that furthermore the magnetic states of the substrate and the substrate-adsorbate system must be carefully controlled in order to obtain physically meaningful results.
Due to the use of the broken-symmetry approximation, most wavefunctions suffer from spin contamination. They are not eigenfunctions of the total spin operator, but contain contributions from other spin states. An energy correction with respect to spin contamination, as discussed by Tada et al., [19] is not available in the CRYSTAL17 program. However, the spin contamination can be calculated as the difference of the actual S 2 -value yielded by the calculation and the ideal S 2 -value of the system. The results are shown in Table 2. The spin contamination decreases with increased ferromagnetic ordering of the spins. The AFM states show large spin contamination of about 8.5-9 in both phases. This is a result of the single-determinant ansatz to model multiconfiguration electronic states. Nevertheless, the addition of the keyword FDOCCUPY results in a  (18) is in all cases the most stable state and has the smallest spin contamination. Therefore, this approach was used to investigate various positions of the adsorbate CO to find the most stable adsorption state. The initial positions on the surface V-and Moatom are shown in Figures 3-6. They are equivalent on both polymorph surfaces. Similar to a previous study of CO and CO 2 adsorption on doped TiO 2 rutile, [30] we tested vertical adsorption C-down on V and Mo, side-on adsorption (side), a tilted structure (tilted) and a rotated structure (rotated). The resulting adsorption energies as well as the optimized distance between the adsorbate and the surface d ads are shown in Table 3. In Figures 3-6 the initial adsorption structures are shown for side, tilted and rotated positions. As in our previous study, [13] the slabs are symmetric in order to avoid artificial dipole moments. Therefore two CO molecules are adsorbed on the top and bottom layers.
All initial CO positions C on V/Mo tilted, side and rotated converged to the same, essentially vertical, adsorption structure shown in Figure 7. The O atom of CO only weakly interacts with the metal atoms of the surfaces, Table 3. Generally, the CÀ M distance d ads is larger on Mo than on V. d ads is about 2.3-2.5 Å on V and about 2.6-2.7 Å on Mo. On the R surface, the adsorption on Mo is the most stable adsorption position with E ads = À 60 kJ/mol. The adsorption on V and Mo on the M 1 surface yield similar adsorption energies of À 49 and À 43 kJ/ mol, respectively. Projected densities of states (pDOS) were calculated for all of the calculated adsorption structures (see Figure S1-S10 in the supplementary material). Only the α (spinup) frontier orbitals are shown which are occupied by the unpaired electrons of V and Mo. The calculated pDOS of the surfaces are similar to the pDOS obtained with sc-PBE0 in a  previous study. [13] The fundamental band gap calculated with M06-D3 is larger compared to the sc-PBE0 results, due to the higher amount of Fock-exchange (27 % compared to 12.7 %). CO adsorption leads to small changes of the surface pDOS, as expected for physisorption. The Fermi energy is slightly shifted to higher values. Small contributions of the C-orbitals are found in the frontier orbitals. The R phase shows a slight localization of the V-d-orbitals near the Fermi level for adsorption of CO on the V-atom. This is not observed in the M 1 phase, where the Vd-orbitals near the Fermi level are already localized. The increase of the band width of the metal d-orbitals due to CO adsorption is more pronounced when the surface is calculated using fixed spin(18) w/o FDO (see Figure S3). This effect is probably responsible for the observed scattering of the calculated adsorption energies. Up to this point, the adsorption energies were calculated in the usual way by subtracting the energies of the separated systems from the energy of the adsorbate system (Eq. 1).
The most stable C-down vertical adsorption position was calculated with a different spin control ansatz to obtain more accurate adsorption energies, because the atomic spin densities still differ between the bare and CO-covered slab, especially in the R phase (see Table 5). We started from the minimum adsorbate structure shown in Figure 7 and gradually increased the CÀ M distance, keeping the vertical CO orientation fixed. For every point the CÀ O distance was relaxed, but the atoms of the slab were fixed at the initial minimum positions. Test calculations showed that relaxation of the surface atoms led to changes in the spin configuration, in particular for the R phase. Applying these constraints it was possible to stabilize the atomic spin distribution over the full distance range between the minimum structure and the separated systems with CÀ M = 7 Å. The adsorption energies in this second ansatz are defined as difference between the energy plateau at CÀ M = 7 Å, and the minimum energy. This is much more consistent in terms of spin configuration conservation than the usual reference of separated slab and CO according to Eq. 1. The resulting potential curves (not corrected for BSSE) are shown in Figure 8. The counterpoise correction was calculated only for the minimum structures of the potential energy curves (see Table 4). It can be safely assumed that the separated systems have no BSSE. For comparison, the previous adsorption energies of Table 3 (now denoted as E prev ads ) are also shown. The CO adsorption energies obtained from the potential curves of the M 1 surface are similar to the previous results. After counterpoise correction similar adsorption energies of À 56 and À 52 kJ/mol are obtained above V and Mo, respectively. The adsorption energies are slightly lower than the previous results of À 49 kJ/mol on V and À 43 kJ/mol on Mo. For the R phase, an adsorption energy of À 46 kJ/mol is calculated for the adsorption on the V-atom. This result is significantly lower than the À 23 kJ/mol yielded previously. The adsorption energy on the Mo-atom changed from À 60 kJ/mol to À 43 kJ/mol. These differences can be explained by the spin distribution changes in the bare and CO-covered slab when using Eq. 1. The much higher consistency of the spin configurations obtained with the  Table 4. Adsorption energies E ads of CO on V and Mo atoms of the R (110) and M 1 (011) surfaces in kJ/mol obtained as difference between minimum and plateau energies in Figure 8 and after counterpoise correction of the minimum structures; M06-D3 results.

Phase
E ads E prev ads on V on Mo on V on Mo M 1 À 56 À 52 À 49 À 43 R À 46 À 43 À 23 À 60 Table 5. Mulliken spin population of the V and Mo atoms of the 4-layer R (110) and M 1 (011) slab without (Surf) and with adsorbed CO (Surf + CO); M06-D3 results obtained with FDO + fixed spin (18). The atom, on which CO is adsorbed, has been highlighted.
Rut 1st layer 2nd layer second ansatz is seen in Tables 6 and 7. Changes of the atomic magnetizations between minimum distance and quasi-separated systems are negligible.
With the second ansatz, we found a difference of � 10 kJ/ mol between the adsorption energies on the R and M 1 phases (see Table 4). The equilibrium CO partial pressure on the surfaces was estimated to investigate the possibility of exploiting the phase transition for CO removal. This was done using the equation based on the Gibbs-Helmholtz: [31] where G ads was approximated by E ads . The CO partial pressure p was calculated at various temperatures for the adsorption on V. In real systems the influence of the dopant Mo, which is present only in small amounts, is assumed to be negligible. The estimated partial pressures in bar are shown in Table 8. At 273.15 K, the partial pressure needed to desorb CO from the R surface is two orders of magnitude higher than on the M 1 surface. One order of magnitude difference between the phases    is still present at 473.15 K. This is an indication that CO removal after M 1 !R phase transition is possible. A reaction with CO as the product may take place on the M 1 surface. Afterwards, the temperature can be raised so that the phase transition to the R phase takes place, where CO could more easily desorb. It has to be noted, that these properties could change with increased CO coverage. Further research on this issue is required. As mentioned above, the accuracy of slab calculations is limited due to the use of a single-determinant wavefunction. With the proposed approach of applying restrictions to the slab geometry and spin configuration with FDO and SPINLOCK in sequential calculations of a potential curve, the corresponding errors cancel out to a large extent. But still the reliability of this approach needs to be assessed by comparison with higherquality methods. Therefore, the CO adsorption was calculated with embedded cluster models at NEVPT2 level. The resulting adsorption energies, again obtained from potential curve scans, are shown in Table 9.
The NEVPT2 results show a strong dependence on the basis set quality. This is mainly due to the basis set incompleteness, the NEVPT2 results therefore include a BSSE, which could not be corrected. The results with the larger def2-QZVP basis set where BSSE is usually small, are more considered as more reliable. The relative stability of CO on R and M 1 obtained with NEVPT2/def2-QZVP and M06-D3 is similar, although the absolute values are larger by � 30 kJ/mol with the former.
The NEVPT2 results give insight into the orbitals that participate in the V-CO interaction. The Mulliken population analysis showed that mainly the V d x 2 À y 2 orbital is occupied in both phases. This is in accordance with the results of the surface models with M06-D3. Overall, the comparison with NEVPT2/def2-QZVP results confirms the validity of our second ansatz.

Conclusions
We investigated the adsorption of CO on the most stable surfaces of the R and M 1 phases of VO 2 . It was found that controlling the atomic magnetizations in order to identify the most stable spin configuration is crucial to obtain numerically stable and therefore reliable adsorption energies, within the limitations of DFT. For both surfaces a ferromagnetic state was found to be the most stable spin state, which also has the smallest spin contamination within the single-determinant approximation. Control of the initial d orbital occupation was found to further reduce spin contamination. Restricted-geometry CO-surface potential curves were calculated for the most stable spin state. By fixing the slab geometry, its spin state could be stabilized which led to numerically stable adsorption energies, that were semi-quantitatively reproduced with a multireference approach. We suggest this procedure for the study of physisorption on open-shell systems, which is particularly suited for atom-centered basis sets. The most stable CO adsorption position on both phases was vertical C-down on the transition metals. An adsorption energy of À 46 kJ/mol on the R phase and À 56 kJ/mol on the M 1 phase were obtained on V. A similar difference is calculated with NEVPT2 and quadruple-zeta basis sets. Thus it may be possible to exploit the R $ M 1 phase transition in a catalytic reaction where CO is a product or an